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Abstract 

This paper develops and analyzes a generic method for reconstructing so- 
lutions to the abstract Cauchy problem in a general Hilbert space, from noisy 
measured data. The method is based on the relationship between a partial 
differential equation and its adjoint equation with control. We demonsrate 
the capability of the method through analysis and numerical experiments. 

1 Problem Description 

Let X be a Hilbert space endowed with the inner product {■,-)x = \\ ■ \\x and 
let A : D{A) C X X he the infinitesimal generator of a strongly continuous 
semigroup 5t : M — > C{X). We are concerned with the following problem. Given 
a C'^ function / : [0,t/] — > X determine the initial condition, xq, of the Gauchy 
problem 

-[t)^Ax{t)+f{t) (1) 

satisfying the measurement condition 

y{t) = Cxii). (2) 

The bounded operator C : X Y retains information about the solution, which 
may only be a portion of the solution, such as boundary values. It is desired to 
determine the initial condition xq given the incomplete (and possibly noisy) mea- 
surements y(i), < t < tf. The method developed in this paper is capable of 
forecasting future states as well, however, we focus on the inverse problem of de- 
termining x(0) due to its practicality and the necessity of dealing with this more 
difficult problem. We are especially interested in the case of partial measurements 



(i.e., the measurements are sparsely distributed over the domain ft). In the fol- 
lowing section, we develop methods for determining the initial state x{0) and we 
demonstrate how the same methods are applicable to forecasting the state x{tf) 
with only minor adaptations. We assume the Hilbert space X is separable, so that 
there exists a complete orthonormal sequence {'Pk}'kLo in X. The approximation of 
x(0) = a;o is given by the truncated (generalized) Fourier series 



where the coefficients satisfy ak — (xq, (pk)- Thus, within this framework, the prob- 
lem reduces to estimating the generalized Fourier coefficients of xq. The problem 
of identifying the initial condition of the abstract Cauchy problem has been widely 
studied. Methods concerning this problem have been covered by Auroux and Blum 
[1], Ito et al [10], and references therein. The monograph by Isakov [5] covers inverse 
problems for PDF in detail. As with many inverse problems, there is extreme dif- 
ficulty in recovering a function from partial and noisy measurements, thus suitable 
regularization is necessary. The main focus of this paper is on the reconstruction 
method coupled with the multi-parameter Tikhonov regularization [6,7]. 

2 An Adjoint Method for approximating the Fourier 
expansion of the initial condition 

In this section, we develop and analyze a new approach for estimating the initial 
condition of the abstract Cauchy problem (1) from time-series data. The method 
developed here involves an indirect computation of the generalized Fourier coeffi- 
cients, based on the adjoint equation of the Cauchy problem. This method has a 
direct link with optimal control theory. Given noisy data, the accuracy and stability 
of the method will be demonstrated in this paper. The general framework of our 
method allows any PDF formulated under the linear semigroup theory to fit into 
this framework. Not only that, but it will be shown that the method can also be 
applied to the less ill-posed problem of forecasting future states of the system. Thus, 
our method may be especially beneficial for applications such as weather forecasting 
or financial futures, where it may be necessary to go both backward and forward. 
Consider the adjoint equation of (1), given by 



where C* G C{Y, X) corresponds to the adjoint of the observation operator C, and, 
likewise. A* is the adjoint of the infinitesimal generator A. Here, u G L'^{Q,tf\Y) 
denotes a control or input to the system. It will be demonstrated that a suitable 
control can be determined for which the generalized Fourier coefficients can be 
approximated by a combination of the control, m, and the data, y. 
Recall the state equation is given by 



m 




^(t)^A*p{t) + C*u{t) 



(3) 



[t)^Ax{t)+f{t) 



(4) 



and the measurements, satisfying 



y{t) = Cx{t), 0<t<tf 



(5) 
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are given, for a known source /. Multiplying (4) by p, (3) by x, subtracting and 
integrating over (0,^/) yields 

j^{x{t),p{t))dt = I {{Ax,p) - {A*p,x) - (C*u,x) + {fit),pm dt (6) 
which implies that 







{x{tf),p{tf))x - {xiO),p{0))x = J {f{t).p{i))x - {u{t),Cx(t))Y dt (7) 



yielding the relation 

*/ 

{x{tf ),p{tf))x - (x(o),p(o))x = j {u{t),m-y{t))Ydt, (8) 



where 



f CSt-sf{s)ds. 
Jo 



The relationship (8) forms the foundation for approximating {xQ,(pk)- 

We recall that the unique mild solutions of the abstract Cauchy problem and 
its dual, with conditions a;(0) — xo,p{tf) — pt^, are respectively given by 

t 

xit) = Stxo + J St-sfis) ds; p{t) = Sl^._^pt, + j S:_^C*u{s) ds (9) 

for each t G [0,t/]. For reconstructing the initial state xq, we assume the controlla- 
bility of the adjoint system (3), which is equivalent to the observability of (l)-(2) 
(i.e., the pair (A, C) is observable). The pair {A, C) is observable if for all a; G X 

*/ 

\\CStx\\ldt>-i\\x\\'' (10) 





for some 7 > 0. Note that by the observability assumption (10), the equation 

y — Mx 

admits a unique solution for y G R{M.), where M. is the operator defined by 

M:^CSt 0<t<tf. 

Furthermore, this unique solution depends continuously on y. The details of infinite 
dimensional control theory are covered in [3]. Having assumed the controllability 
of (3), we define the operator ^ : L'^ {0, tf;Y) ^ X hy 



-1 

:= J S*C*u{s)ds 



(11) 
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for u G L^{0,tf;Y). By construction, the adjoint equation (3) evolves backwards in 
time. If p{tf) = 0, then the adjoint satisfies 



p(0) = / S*C*u{s)ds. 
Jo 

By the controUabihty/observability assumption, we know a unique solution to 

^u^piO), (12) 

exists for p{0) £ R{^). However, in practice, the exact controllability of (3) is, in 
general, not true, so we assume the condition (11) holds approximately, i.e., there 
exists Ue such that 

\\^u,-pm\x<e (13) 

for any e > 0. Whenever this relationship does not hold (or holds only approxi- 
mately), we must suitably regularize the problem, so that a reasonable u can be 
obtained. We are interested in solving for w, since by relationship (8), we can obtain 
the /c*^ Fourier coefficient of xq as 

{xo,^k)x = / {u{t),y{t) - ^{t))Y dt, 
Jo 

whenever pt^ =0. 

We proceed by defining a collection of adjoint functions pk (0) = (pk , such that 
{V'felfcLo forms an orthonormal basis for a finite-dimensional subspace Xm C X. 
Then {(a:;(0), <(5fe)}feLo ^''^ generalized Fourier coefficients for x{0). By the con- 
trollability assumption (10) and by utilizing relation (8), we can determine the 
Fourier coefficients of a;(0) by solving the operator equations 

^Uk ^Vk, 0<k<m. (14) 

If (11) or (13) holds, we will construct stable approximations Uk using a suitable 
regularization method. An example of such a regularization method for determining 
one-dimensional Uk is to solve the minimization problem 

min llifu-^lll+ryi f ' \u{t)\ dt + \u'{t)\' dt, 

where the first term corresponds to the sparsity of the approximate solution Uk(t), t G 
[0,tf], while the second term corresponds to the smoothness of Uk- We note that the 
smoothness of Uk may affect noise dampening (see Remark 2). Such regularization 
methods are described in detail in the papers [6], along with criteria for selecting 
the regularization parameters ?7i , 772 • 

Our approach is based on the fact that for each basis function ipk there exists 
a control Uk G i^(0,t/;F) such that ^Uk = ^k (or W^u^ - p{0)\\x < £)■ The 
controls Uk(t) are determined in such a way that each adjoint pk is driven from zero 
at time tf to Pk(0) = Vk, for a suitably chosen (pk- With each Uk determined, we 
construct the approximation for xq by 

^o" ^Yl ^^"^^ 

k=0 
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where the generahzed Fourier coefhcients are approximated by 

{xo,ipk)x ^ 1 {uk{t),y{t) - ^{t))Y dt ^ Uk (15) 
Jo 

using the relation (8) and equation (14). 



Further analysis of the method is detailed below, including the error analysis in 
Theorem 2.1. The following summarizes the method for estimating xq. 



Dual Method for reconstruction of xq: 


1. Pick an 


orthonormal basis, {fkYk-o -^"^ ^ 


2. For each k solve ^Uk = (fik to find Uk G L^{0,tf; Y) 


3. Form the estimate for xq, 




m 
fc=0 


where 


ak= / {uk{t),y{t)-^it))Ydt 
Jo 


with 


at)= f cst-sf{s)ds. 

Jo 



The well-posedness of the method follows from the controllability assumption 

r \\CStxfdt>-i\\x\\\ (16) 
Jo 



for all X £ X. 

Using the method for forecasting a future state 



Now, we briefly introduce how the method is utilized for the purpose of forecasting 
a future state x{tj). For this purpose, we assume the adjoint (3) is null-controllable, 
i.e. there exists u G L^(0, t/; F) such that p(0) = and 

^u = -St,p{tf). (17) 

Recall that p evolves backwards in time (with respect to the evolution of x) . In gen- 
eral, the exact null-controllability may not hold, however we assume the condition 
(17) holds approximately, i.e., there exists such that 

\\^u, + St,p{tf)\\x <e, 

for any e > 0. With Uk determined, the generalized Fourier coefficients are approx- 
imated by 

{x{tf), ifk)x = / ' {ukit), m - y{t))Y dt 
Jo 
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where Uk is the approximate sohition to 

ifufc = -Stf(pk- 

For the final state case, the method is well-posed under the assumption of null- 
controllability of the adjoint control system, i.e. 



The method is summarized as follows 



Dual Method for reconstruction of xtfi 

1. Pick an orthonornial basis, {v'felfeLo fo^' -^m C X 

2. For each k solve = —Stffk to find Uk G L^{0,t f;Y) 

3. Form the estimate for xt^ , 

m 

m \ ^ 

Xtf = 2^akfk 

k=0 



where 



with 



ak^ / {Ukit),m -yit))Ydt 
Jo 

m= f cst-sf{s)ds. 

Jo 



'i'he novelty of this method is, in part, due to the fact that it is not necessary 
to compute the time history of the adjoint, p. However, the method utilizes the 
information available from the adjoint in order to accurately reconstruct xq. By 
utilizing the norm, we are able to construct sparsely distributed controls, which 
can aid computational efficiency. Furthermore, the method is quite robust to noise, 
as the actual inverse problem does not involve the noisy data. 

Remark 1 We also note that there is a stochastic interpretation of this method. 
Assume x,p are random variables satisfying the linear stochastic differential equa- 
tions 

dx = {Ax{t) + f{t))dt + adBu -dp = {A*p{t) + C*u[t))dt (18) 

where Bt is the Brownian motion, and a is the standard deviation ( diffusion coef- 
ficient). Then, by the relation (8) we have 

{xo,ipk)x= [ \uk{t),y{t))Ydt- f \f{t),pk{t))xdt + a f ' pk{t)dBt 
Jo Jo Jo 

which implies that 

E[\{xo,ipk)x - r {uk{t),y{t)-m)Ydt\']^E[a^\ f pk{t)dt\% 
Jo Jo 

Thus, the mean square error in approximating the Fourier coefficients is propor- 
tional to the standard deviation, a, of the Brownian motion, regardless of that fact 
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that p(tf) — (in the case of estimating xq). Determining the control, Uk, can be 
cast as 

min ll^u - +/?crM \p{t)\'^ dt < k < m 

u&L^O.tfX) Jq 

where 

p{t) = j'^' SUC*u{s)ds. 

Thus, we select the parameter j3 so that e^ + /3(T^ is balanced, where e is the accuracy 
of the fidelity term 

^Uk - ipk = £■ 

The following theorem provides the error estimate of our reconstruction method 
in the real Hilbert space setting, as well as justification for the method based on 
mixed regularization. In short, there are two sources of error in approximating the 
Fourier coefficients. The first source of error is due to the ill-posedness of J^Uk — fk, 
while the second source of error is due to the noise, S, in the observed data. The 
errors must be balanced to obtain the best possible solution. The proof is omitted, 
as it is a straightforward application of the Cauchy-Schwarz inequality. 

Theorem 2.1 (Error Estimate) Suppose {A*,C*) is approximately controllable, 
there exists Uk G L^{0,tf]Y) such that 

W^Uk - >Pk\\x < £k 

for each < k < m. If we define, 

v{t) = y\t) - y{t) 

and 
then 

m 

\\xo - xTWx < lk-0 - a;™IU + E (e^H^olU + c(<5, i/)||wfc(t)|U) 

k=0 

where Z = L^{0,tf;Y) and 

Wxo-x'^Wx 

is the truncation error of the generalized Fourier series. Furthermore, if xq G 
C'in), then 

\\xo - xTWx < E i^Mx + c{S,tf)\\ukmz) (19) 

fc=0 

for m sufficiently large. 

Better error estimates may be realized, however, the results of Theorem 2.1 also 
provide justification for the regularization methods. By the estimate, 

m 

\\xo - xTWx < \\xo - x"'\\x + II {{^Uk - ipk,xo)x + {uk{t),v{t))z) ^fikllx 

k=0 
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we immediately see the need for appropriately solving Uk- If the noise level, S, is 
large we must obtain controls which arc sufhcicntly regular, so that the term 

{uk{t),v{t))z 

is small, while simultaneously ensuring \\^Uk — ^k\\x is small. The following remark 
further justifies the previous statement. 

Remark 2 Suppose the noise in the data is highly oscillatory, such as cos{l'Kt). 
Then the error in the Fourier coejficients has the term 

1 1 

Uk{t) cos{lnt) dt = — [ u'f.{t) sin(/7rt) dt. (20) 
In J 



That is, the highly oscillatory parts may be damped by In, if Uk is sufficiently smooth. 
Thus, we utilize a penalty which enforces smoothness on the control Uk- 

It is also apparent that the accuracy, e/c, in solving 

is necessary for an accurate reconstruction of xq. In practice, we must balance the 
accuracy of solving ^Uk — (fik and the regularity imposed on Uk via the regular- 
ization methods. This concern is addressed in Section 3.1 where we discuss how to 
balance the method to obtain stable but accurate solutions. 



2.1 Variation of the Dual Control Method 



In this section, we outline an alternate procedure for obtaining reconstructions of 
the initial condition, xq. This approach is based on the adjoint control approach 
developed in the previous section. Rather than selecting a collection {pfc(0)}^g to 
be a basis for X, we select {uk{t)}^^Q to be a basis(not necessarily orthonormal) 
for Z = L^{0,tf; Y). Assuming the relation (11) holds, we construct the adjoint set 
{PfclfcLo by the relations 

^Uk = Pk- 

Note that the collection {pk}k=Q linearly independent under the assumption that 
{A,C) is controllable, i.e., 

R{^) ^ N{^) = 0. 

Thus, if {A, C) is exactly controllable, we form an orthogonal(orthonormal) basis 
by the Gram-Schmidt method. The coefficients of xq are computed by defining the 
Gram matrix 

Gk,i = {pk,Pi)x 
and setting f3 — (/3o, . . . , /3m)* such that 



/ 



13 ^G- 







\ 



y')dt 



dt. 



\ 
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The coefficients f3 can be computed efficiently by the Cholesky decomposition G = 
LL*, since G is symmetric positive definite. Again, the algorithm is well-posed 
under the exact controllability (16) of the adjoint system which, in general, may 
not be true. If the adjoint system is not exactly controllable, care must be exercised 
to ensure the set {^Uk}^^Q is linearly independent. 



Variation of Dual Control Algorithm: 

1. Pick a basis {ufc(<)}^Lo ^r [/„ C L'^{0,tf;Y) 

2. Compute pk by Jfufe = pk 

3. Compute the Gram matrix Gk.i = {pk,Pi)x 

tf 

4. Set Hi. = J {uk{t) , S,(t) — y{t))Y dt and compute the approximate 



Fourier coefficients (3 = G~^y 

5. Compute the approximation 

m 

x^ = Y.Pk^Uk 

k=0 



There are several potential advantages to this approach. Namely, one can directly 
regulate the properties of the controls Uk , such as smoothness or sparsity. Secondly, 
the operator ^ does not need to be inverted. However, since the pair {A*,C*) is 
not necessarily controllable, we are not guaranteed linear independence of the set 
{Pk}T=o- Thus, solving 

/ J^'{uo,y'}Ydt \ 
G/3 - ; (21) 

\ Io'{U'm,y^)Ydt J 

for (3 requires regularization. This method only requires the solution of one ill- 
posed problem, but requires the formation of the m + 1 adjoints pk- Therefore, 
this method may be less expensive than the dual control method, however, with a 
tradeoff in accuracy. 

2.2 Implementation Issues 

In this section, we discuss the necessary numerical issues for the implementation 
of the methods developed in this section. For the numerical implementation for 
solving the dual control problem we use the Crank-Nicholson scheme 

-^-—^=A*i^ — :P^ + C'*Uk+i (22) 
At 2 ''+2 ^ ^ 

for (3) where Uk+1/2 is evaluated at the mid-point of the interval [tk,tk+i] and 
tk — ktfAt. At the time step fc -I- 1 the solution is computed by 

/+i = (/ + ^aA " (/ - ^aA / -At(l+ ^aA " G*u,^. (23) 
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utilizing the (1, 1) Pade approximation for exp{—A*At). In the dual control formu- 
lation, the discretized problem for each control u is formulated as 

min WL^u (p\\y + l^rpiu) 

where L" — (M")*, given 

^ [C,Cn.i{AAt), . . . ,Cn.i{AAtr-^] (24) 

and Um is a finite-dimensional subspace of L^{0,tf;Y) and ijj is a chosen penalty 
term. 

If necessary, higher order Pade approximations may be considered, which are of 
the form 

Pm , , ao + aiZ + . . . + flmZ™ 

= = — ■ , „ (25) 

Qn bo + biZ + . . . + bnZ" 

where the degree of P, Q is not more than m, n respsectively. Higher order Pade 
approximations of semigroups are discussed in detail in the paper [14]. 

Operator Splitting for Convection-Diffusion Equation 



The Crank-Nicholson scheme works well for the diffusion dominant case, however, 
for the convection dominant case it is necessary to solve the problem more accurately 
( and such that the physics are obeyed). In this section, we describe the numerics 
for the initial condition estimation of the convection-diffusion equation 

^{x,t)=c{x)-Vv{x,t)+V-{d{x)Vv){x,t)+f{x); v{x,0) = vo{x) (26) 

where c{x),d{x) are the convection and diffusion coefficients, respectively. 

The reconstruction methods have a natural extension to such problems, using a 
differential operator splitting 

^^^Lv{t)^{A + B)v{t) 

where A,B e C{X) 

For the numerical solution of the convection-diffusion equation, we consider two 
stage Strang operator splitting 

v{x,t + At) = S'LSl^S'Lv{x,t) 

2 2 

where Sf , are the semigroups corresponding to the parabolic and hyperbolic sub- 
problems, respectively. That is, Sf, are the Co-semigroup semigroups generated 
by A, B respectively. 

Assuming a constant convection coefficient c, we solve the hyperbolic subprob- 
lem via the method of characteristics v{tn+i,x) = v{tn,x — cAt) where the right 
hand side is evaluated via cubic interpolation. As in the previous section, we use 
the Crank-Nicholson method for solving the parabolic subproblem, using the ap- 
proximating polynomial 

2 + z 

for the approximation of exp(v4Ai). 
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3 Generalized Multi-parameter Regularization for 
Control Solution 



This section is devoted to discussing the multi-term regularization method utilized 
for solving (11), without going into detail. In general, rather than solving 

directly, we seek a minimum of 

Jp[u) = <j){u,^)+p^{u), (27) 

over u G C, where the fidelity term, 0, is chosen based on the noise statistic, 
while ip is chosen based on which class the solution x should belong to. Whenever 
0(m, ip) ~ \\J^u — ip\W and tpiu) = this formulation coincides with the classical 

Tikhonov regularization 

mm\\\^u~^\\\ + l-\\u\\\. 

The main drawback to this method is the single regularization term ip. Modern 
day scientific problems typically involve applications where the standard Tikhonov 
regularization fails to capture the full set of distinct features in the physical solution. 
Many research efforts have been devoted to improving the standard regularization 
techniques for a wide range of applications (see [2,8,16,17] for example). It is not 
a goal of this paper to cover this in detail. These references and the references 
therein provide a thorough study of such methods. Especially in the field of image 
processing, the solution often exhibits a multiscale structure typically described by 
multi-resolution analysis. In such applications, single parameter regularization can 
oversmooth the solution such as the case of ■0 = || ■ 11^2 or exhibit stair-case effects 
such as the case of -0 = || • Htv- In order to capture the multiscale structure of 
solutions without introducing oversmoothing or staircasing, many research efforts 
have focused on mixed regularization approaches, such as combining the penalty 
term with the TV penalty: 

min \j\^u~ (^p d^+'^J d^ + mj | Vu| (28) 
a n n 

To capture multi-scale solution profiles, we employ the multi-parameter Tikhonov 
regularization technique, i.e., we minimize 

Jr,{u)^<j){u,ip)+rj-xP{u). (29) 

The terms 0, xp are known as the fidelity and regularization terms, respectively. 
Here, {ipk}k=i is the set of regularization terms, {r]k}^=i £^re the regularization 
parameters, and we take the dot product 

n 

V ■ i>iu) = ^?7fe0fe(u) 

k=l 

for ry = (r/i,7]2, . . . ,?7„) and ■0(u) = ('0i(u), 02(w), ■ • ■ ,0'n(w))- The functionals 4>,ip 
can be chosen based on any a priori information about the problem and its exact 
solution. Then, 

Urf = arg min J"^ (u) 
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is taken as the regularized solution. For instance, in the case of a multiscale image 
with a smooth region and a stepped region, one may consider the L^-TV regular- 
ization (28). In this work, we also consider the penalty term 

,Piu)^\\u\\l p<l, 

to enforce sparsity in the solution. 



3.1 Balance principle 

We discuss here the balance principle for the single-term regularization, based on 
the paper [6,11]. Prior to this selection rule, most selection rules (e.g. Morozov's 
discrepancy principle) were based on either the performance level (noise) 

(/)(U,V3) 

or the complexity level 

alone. The selection rule developed in [6, 11] is based on balancing the perfor- 
mance level and the complexity level. Consider maximizing the conditional density 
p{{u,T, X)\(p) ~ p{ip\{u,T, X))p{u,T, \) where (r, A) are density functions for (f>,'ip, 
respectively, both having Gamma distribution. The balancing principle is derived 
from the Bayesian inference [11] 

min T(f){u,(p) + XtP{u) + fjoX — dolnX + rjiT — aihiT. (30) 

{U,T,\) 

By definition, (u. A*, r*) e X x M+ x M+ is a critical point of (30) if 
u* = arg min{0(u, ip) + X* {T*) — l^p{u)} 

u 

(j){u*,ip) + m ~ ai:^ ^0 

By optimality, the regularization parameter satisfies 

, ao(j){u*,ip)+rii 
^ = 7T^^ ■ ^^^> 



The authors arrive at the selection criterion 



< d < 1, 



by rescaling ao as ctq in order to ensure the conditions 

2 

lim 77(0-0) = 0, lim = 

o-o^O o-o-i-0 77(0-0) 

are satisfied, where CTq the variance. Further discussion on the validity of this 
method, as well as the selection of the constants ai,rio,d, can be found in [6]. The 
following iterative algorithm for determining 77 is utilized: 
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Iterative algorithm to solve for {ur^,ri): 

Choose an initial guess /3o > 0, and set k = Q. Find (u^, Pk) for fc > 1 as 
follows: 

1. Solve for Uk+i by the Tikhonov regularization method to obtain 
Uk+i = arg niin = {(j){u, ip) + I3kip{u)}. 



2. Update the regularization parameter by 

Pfc+i = a—, r- ■ 

^(uk+i) + -no 

3. If a stopping criterion is met, stop, else set A; = fc + 1 and repeat 
from step 1. 



4 Numerical Tests 

4.1 1-D Diffusion Equation 

In this section, we consider inverse problems involving the 1-D diffusion equation 

v{Q,t)=Q^v{l,t) (32) 



y{t) = Cv{t) 

with Dirichlet boundary conditions, where the measurements are restricted to a 
subinterval Vlg dVl — [0, 1], for the time < i < 1. Specifically, the operator C takes 
average measurements over the two intervals (.23, .31) and (.46, .53). The thermal 
conductivity, d, is potentially variable in space, but known. The 1-D diffusion 
equation is formulated as an abstract Cauchy problem (1) where 

dx dx 

and 

dom(A) — { u G i^(0, l)|t;, ^ are absolutely continuous, 
^ ei2(o,l) and«(0) -0-t;(l)}. 

It is a standard exercise to show that A generates a strongly continuous semigroup. 



Example 1 : Spatially varying diffusion coefficient 

For this example, we consider the case when the thermal conductivity is spatially 
variable. In particular, we take 

d{x) = 1.0625 - (x - i)4 
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and the initial condition is given by 

z;o(x)-e-20o(.-i)^^ 

We take the basis ipk — {sin(fc7rx)})!LQ and solve for the controls, Uk, using the 
L^-H^ regularization method, for m = 8. It should be pointed out that the ab- 
stract Cauchy based dual control method does not make any assumptions on the 
coefficients of the PDE. Assuming a noise level of 10%, we obtain an accurate and 
stable reconstruction with the parameters 771 = 5 x 10~*, ry2 = 1 x 10""'^°. The 
corresponding results are depicted in Figures 1. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(a) 



- - Appraximale solution I 
Exact solution 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

(b) 

Figure 1: (a) Thermal conductivity; (b) Reconstruction with 10% noise in measure- 
ments versus Exact initial condition. 



Comparison of basis choices 

Here, we compare the reconstructions obtained by two different basis choices. For 
this example, we take 

5(a;-i)4, 0<a;<i 
1 1 <r T < 1 



d{x) 



1^ 

16 
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as depicted in Figure 2, and the initial condition is given by 

^;o(x)=e-200(.-i)^^ 

We solve for the controls Uk using the L^-H^ regularization method. Assuming a 
relative noise level of 5%, we obtain an accurate and stable reconstruction with the 
parameters rji = 5 x 10~'^,r]2 = 1 x 10~^^, by computing m = 8 coefficients, using 
Daubechies wavelets. In Figure 3, one can see a comparison of two reconstruc- 
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Figure 2: Thermal conductivity. 

tions using a standard sine basis and Daubechies-10 wavelets. The reconstruction 
obtained using the Daubechies wavelets is more accurate and stable than the sine 
basis reconstruction, even with well-tuned regularization parameters. This example 
illustrates how the basis choice affects the resulting reconstruction. 

4.2 2-D Diffusion Equation 

In this section, we consider inverse problems involving the 2-D diffusion equation 

^ix,t) =dAvix,t) + f{x) 

u(a;,t)=0 xedn 

v{x,0) = vo{x) 

where x — {x,y) E C M^^ . As in the 1-D case, we work on the time interval 
< i < 1. 

The 2-D diffusion equation can be cast in the abstract Cauchy framework where 
A coincides with the closure of the Laplace operator, defined by 

for every / in the Schwartz space 

^(M") := (/ e C°°(R") : lim = for ah A; e N and a £ N" j . 
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- - Approximate solution I 
Exact soiution 
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Figure 3: (a) Reconstructed initial condition using Daubechies-10 wavelets with 
771 = 5 X 10~^,772 = 1 X 10^^^; (b) Reconstructed initial condition using sine basis 
with 7/1 = 5 X 10^6, 772 = 1 X 10-^°. 



4.3 2-D Convection-Diffusion Equation 

In this section, we present severeal numerical results for inverse problems involving 
the convection-diffusion equation 

v{x, 0) = vo{x). 

For the results presented here, we assume c{x) = c, d{x) = d are constant (or at 
least locally constant), and we take f{x) = 0. For both simulations, the domain is 
taken as the unit square il = [0, 1] x [0, 1]. 

Initial condition reconstruction 

We consider the initial condition reconstruction problem with d = .l,c = (515) 
known, where the observation operator is defined by 
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where /J,(ris) is the volume of the set ilg- That is, we take average measurements 
over a sample set fig C For this simulation, we take nine measurement locations 
equally spaced over the domain, each location of size x jij. The corresponding 
contaminated measurements are depicted in Figure 4b. Using the operator split- 
ting technique outline in Section 2.2, the convection-diffusion equation fits into the 
abstract framework (1). The exact initial condition is 

and we solve the corresponding inverse problem using the Li-Hi regularization, 
with basis functions 

ipk,iix,y) = sin(fc7ra;) sin(;7ry). 

As can be seen by comparing Figures 5a and 5b, the method for reconstructing 
the initial condition performs well with the parameters r/i = .03, 771 = 1 x 10~^. 
Depending on the basis choice, small errors are expected due to the truncation of 
the generalized Fourier series. In this case, we have small oscillations indicative of 
the sinusoidal basis. 




(a) 
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Figure 4: (a) Nine measurement locations depicted in red; (b) Noisy measurements 
used for reconstruction compared with exact measurements. 



17 



(a) 




y 



(b) 

Figure 5: (a) Exact initial condition; (b) Reconstruction with 10% noise in mea- 
surements. 

5 Concluding remarks 

The abstract Cauchy problem provides a unified framework for the analysis of sys- 
tems governed by PDE. The methods developed in this paper allow for the system- 
atic reconstruction of initial conditions of the abstract Cauchy problem. In particu- 
lar, the dual control method coupled with the multi-parameter regularization yields 
a method that is very tunable and robust. By an appropriate basis selection for 
the problem at hand, and by selecting the parameters in the regularization frame- 
work based on the balance principle, a reconstruction filter is determined based on 
the governing PDE. Depending on the problem size, there may be significant over- 
head in computing the controls (11). However, once computed, the controls can 
be banked (or stored) for future use. Thus, if one carefully selects the basis and 
the parameters are tuned to the noise and a priori information about the solution, 
the method can potentially be implemented in real time, simply by integrating the 
controls against the data. 

Diffusion processes and parabolic equations fit particularly well into this frame- 
work, due to the necessity for stabilizing the dynamics backward in time. The 
method accurately reconstructs both initial conditions and point sources of diffu- 
sion processes, and allows the forecasting of future states. Thus, the tool provided 
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is valuable for problems where numerous calculations are required based on sensor 
data, and for problems where integrating forward and backward in time is impor- 
tant. 

Based on the multi-parameter regularization, the methods developed are par- 
ticularly suited for problems involving a locally supported source, such as point 
sources, as well as those with sparsely distributed data. The sparsity optimization 
works well for both identifying initial conditions/sources that are locally supported, 
as well as for selecting the necessary control profile. 

Certain questions still remain and extensions to more difficult problems can 
be realized. Specifically, nonlinear problems can be treated in a similar manner, 
through the development of nonlinear dual control filters. In a forthcoming paper, 
we describe the nonlinear method for equations such as the one-dimensional viscous 
Burger's equation 




the Korteweg-de Vries (KdV) equation 

and its generalizations (e.g. the Novikov-Veselov equation). We are also interested 
in inverse problems regarding the incompressible Navier-Stokes equations 

divt; = y-^=0 (xeM",t>0), 

with initial conditions 

v{x,0)=vo{x) {xeR"). 

The three-dimensional Navier-Stokes equations have important applications, such as 
weather modeling, aircraft design, rheology, etc. Due to the practical need for con- 
sidering three-dimensional Navier-Stokes, efficiency must be addressed. In this case, 
the solution for the controls must be performed efficiently, though, once computed, 
this framework may be ideal for large scale problems since the filter coefficients can 
simply be banked. Thus, future research for this method also involves addressing 
computational efficiency. 
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